Criticality and Scaling Relations in a Sheared Granular Material 
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We investigate a rheological property of a dense granular material under shear. By a numerical 
experiment of the system with constant volume, we find a critical volume fraction at which the shear 
stress and the pressure behave as power-law functions of the shear strain rate. We also present a 
simple scaling argument that determines the power-law exponents. Using these results, we interpret 
a power-law behavior observed in the system under constant pressure. 
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I. INTRODUCTION 

Soft glassy systems such as foams, colloidal suspen- 
sions, emulsionsjpolymers, glasses Q,Q,E|, and granular 
materials have a strongly non- linear response to 
an external perturbation. In such systems, the relation 
between the stress a and the strain rate 7 characterizes 
the system behavior. Although it is known that the re- 
lations are diverse and specific to individual systems, a 
universal law for a certain class of systems may exist. 

In particular, in sheared granular materials under con- 
stant pressure p, one of the authors (Hatano) has found 
a relation [7j 



with 




(1) 



(2) 



by a numerical experiment using the discrete element 
method. Here, a is the maximum diameter of the par- 
ticles (their diameters are uniformly distributed in the 
range [0.7a, a]) and m is the mass of the particles 0. As 
demonstrated in Fig. ^ the exponent <p is not inconsis- 
tent with 1/5 in the range 10 -3 < I < 10 _1 . Surpris- 
ingly, the power-law behavior given in Eq. Q is observed 
in the cases that p/Y ~ 10~ 3 and 1CP 5 , where Y repre- 
sents the Young modulus of the particle. For example, 
one can experimentally obtain the power-law behavior 
under the constant pressure p = 10° — 10~ 2 MPa by using 
polystyrene with Y = 3GPa. Since / = 10~ 3 corresponds 
to the shear rate 10° — 10 1 /sec in this example, the shear 
condition leading to Eq. is experimentally possible. 
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FIG. 1: a/p as a function of I. This result was obtained 
for a model similar to that explained in the text. The main 
differences are as follows: (i) the top boundary in the z di- 
rection is modified so as to maintain a constant pressure and 
(ii) the shear is applied directly from the moving layer at 
the top and the bottom. (See the inset.) The parameter 
values are as follows: L x /a = L y /a = 25, N = 10800, and 
t]/Vkm = 1.0. pa/k = 1.92 x 10" 3 (data 1), 3.75 x 10" 5 (data 

2) , and pa/k = 4.35 x 10~ 3 with L x /a = L y /a = 15 (data 

3) . Furthermore, the square and circle symbols represent the 
constant pressure data obtained from Figs. [5] and [3] where 



pa/k = fl 
symbol) . 



1.25x10 (square symbol) and 5.4 x 10 (circle 



Stimulated by this result, in the present paper, we 
consider the power-law behavior of stress-strain rate re- 
lations in sheared granular materials by investigating a 
model granular system with the Lees-Edwards bound- 
ary conditions. In this idealized system, we demonstrate 
that there is a critical volume fraction at which the shear 
stress and the pressure (normal stress) behave as power- 
law functions of the shear strain rate in the limit 7 — > 0. 
From these power- law behaviors, we derive the scaling 
relation a/p ~ J 1 / 5 in the limit I — ► at the critical 
volume fraction. Note that this critical condition does 
not correspond to a constant pressure. We then present 
a simple interpretation of Eq. Q for the system under 
constant pressure. 



2 



II. MODEL AND NUMERICAL RESULTS 

Here, we describe our computational model. The 
system consists of N spheres of mass m in a three- 
dimensional rectangle box whose lengths are L X1 L y , 
and L z along the x, y, and z directions, respectively. 
In order to realize an average velocity gradient 7 in 
the z direction and average velocity in the x direction, 
we impose the Lees-Edwards boundary conditions 0. 
The particle diameters are 0.7a, 0.8a, 0.9a and a each 
of which is assigned to N/A particles. When the dis- 
tance between two particles is less than the sum of their 
radii, r\ and an interaction force acts on each of 
them. This force comprises an elastic repulsion force 
k(5r — (ji + r 2 )) and the viscous dissipation force r]5v, 
where Sr and Sv represent the relative distance and ve- 
locity difference of the interacting particles, respectively. 
For simplicity, we do not consider the tangential force 
between the interacting particles. We study the spe- 
cific case where L x /a = L y /a = L z /a, N — 1728 and 
77/ V km = 1.0. The control parameters in this system 
are the volume fraction v = ^2 i _ i waf/(6L x LyL z ) with 
the ith particle diameter dj, and the dimensionless shear 
rate T = 7^/m/fc. We then calculate the dimensionless 
shear stress S = aa/k and the dimensionless pressure (in 
the z direction) II = pa/k. See Ref. as the calcu- 
lation method for £ and EL Note that k/a provides an 
approximate value of the Young modulus of particles. 

We express the dependence of £ and II on (T, v) as 
£ = / CT (r, v) and II = f p {T,v), respectively. Figures [21 
and[3]display these functions with respect to T for several 
values of v ^j- These graphs clearly show that there 
exists a critical volume fraction v c at which the power 
law behaviors are observed as follows: 



/ CT (r,o ~ r a , 
/ p fE> c ) ~ r' 3 , 



(3) 
(4) 



in the limit r — ► The values of the exponents 

will be discussed later. Here, it is worthwhile noting 
that similar graphs were obtained in Ref. [ljj with the 
argument on the effect of finite elastic modulus. Indeed, 
these graphs in this reference suggest the existence of 
the critical state, although the power-law behavior was 
not mentioned explicitly. Upon numerical verification, 
we found that the critical volume fraction corresponds 
to the jamming transition point defined as the volume 
fraction beyond which a finite yield stress appears [l2j| . 
In this paper, we do not argue the nature of the jamming 
transition, but focus on the power-law behaviors given in 
Eqs. © and J2J|. Note that a similar critical state was 
obtained for a sheared glassy system p|. 



III. THEORETICAL ARGUMENT 
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FIG. 2: E as a function of F for several values of v. The 
thick solid line represents E oc r 5//7 that is estimated from 
our theoretical argument. Note that the Bagnold scaling [T3 | 
is observed for the case in which v = 0.630 and F < 10 -4 . 
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FIG. 3: II as a function of F for several values of v. The 
thick solid line represents II oc F 4//7 that is estimated from 
our theoretical argument. Furthermore, the lines I = const, 
are drawn to help us to understand the functional form of a/p 
over 



defined as 



JV 



i 2 ), 



(5) 



where Vi denotes the velocity of the i-th particle. Al- 
though T is not a parameter of the system but is deter- 
mined by r and v, it is considered that physical processes 
in granular systems are described in terms of the kine- 
matic temperature |l4j . In particular, the time scale 
of energy dissipation is assumed to be determined as 
(y/T/m/a) -1 . One can verify the validity of this as- 
sumption by investigating the energy balance equation 
in the steady state [l4|: 



The main idea in our theoretical argument is to con- 
sider dimensional analysis with kinematic temperature T 



N T 3 / 2 



L x L y L z a^/m 



£77, 



(6) 
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which is rewritten as 



T 

k^? 



3/2 



nr. 



(7) 



Figure 0] indicates that Eq. J7J is plausible as the 
first theoretical attempt, although a slight deviation is 
observed. Based on this result, hereafter, we assume 
that the time scale of the energy dissipation is given 
by (y/Tj^/a)- 1 . 
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FIG. 4: u{T/(ka 2 )) 3/2 versus EF for two cases in which v 
is close to the critical volume fraction v c . The solid line cor- 
responds to Eq. J7J. Since we find the best fitting to be 



v(T/(ka')) 
0. 
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(Er) ' , there is a slight deviation from Eq. 



Now, we consider a set of dimensionless constants 
using the energy dissipation ratio. Let us recall that 
we have three independent dimensionless parameters 
(u, r, rj/Vkm). We want to use dimensionless parame- 
ters each of which represents the ratio of a time scale 
of a physical process to that of energy dissipation. It is 
then reasonable that the following two time scales are 
important: the inverse of the shear rate 7 -1 and the re- 
laxation time 77/fc of the elastic displacement during the 
inte ract ion between two particles. Thus, instead of T and 
77/Vfcm, we introduce 



a 

6 



1 

ka 



(8) 
(9) 



Then, we wish to determine the functional forms of 
E(6,6) and 11(6,6) m the limit r — > at the critical 
volume fraction. In order to restrict the possible forms 
of the functions, we further assume that a does not de- 
pend on k. Noting that £2 — ► for T — > 0, the function 
becomes 



£^A(6)6, 



(10) 



Here, we consider A(6) on the basis of the theory de- 
termining the behavior of shear stress near the critical 
state for a dense colloidal suspension [l5|. According to 
this theory, the critical behavior is described by an or- 
der parameter equation that is similar to the Ginzburg- 
Landau equation for magnetization under a magnetic 
field. Assuming that this description is valid for the 
present problem, we write 



co(y - v c )A(£i) - aAfo) 3 = & 



(11) 



with numerical constants cq and c\. Note that, in cases 
of magnetic materials, A(6) and 6 correspond to mag- 
netization and a magnetic field, respectively. Because 
the first term vanishes at the critical volume fraction, we 
obtain 



1/3 



(12) 



Following these assumptions, we can calculate the ex- 
ponents a. Concretely, Eq. (fTUf> with Eq. iQ becomes 



Combining Eq. (fP3*)l with Eq. (7J, we derive 

T : 



(13) 



(14) 



The substitution of this into Eq. ljT3")l with Eqs. (jSJ) and 
© yields 



(15) 



Thus, we obtain a — 5/7 in Eq. @, which is consistent 
with the numerical experiment as shown in Fig. [5] In a 
similar manner, we obtain 




n : 



n 



10/7 



p4/7 



(16) 



under the assumption that 11(6,6) — £2- This assump- 
tion implies that II does not depend on 6 because the 
normal stress is not directly influenced by the shear rate. 
As shown in Fig. 01 /? = 4/7 in Eq. J3J is consistent with 
the numerical experiment. 



IV. INTERPRETATION OF EQ. Q 

We next study the power-law behavior observed in the 
system under constant pressure on the basis of the results 
obtained above. First, from Eqs. © and (0J, which are 
valid at the critical volume fraction v c , we derive Eq. 
with 



where -A (6) represents the direct contribution of V that 
is not expressed through the dependence of T on Y. 



(17) 
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Using the values a — 5/7 and /3 — 4/7, we obtain cj) = 
1/5. This value of 4> is consistent with the numerical 
experiment. However, it should be noted that the scaling 
relation is obtained at the critical volume fraction, not 
for systems under constant pressure. 

In order to discuss quantitatively the behavior of the 
system under constant pressure, we denote the volume 
fraction and the shear stress measured on this system as 
v = <7„(T,n) and E = g a (T,H). We then wish to deter- 
mine these functions from Figs. [21and|3J First, the point 
(r, II) in Fig. ^determines the volume fraction uniquely. 
We assume here that this volume fraction is realized in 
the system under constant pressure II with shear rate 
r and that the shear stress in the system is determined 
by using the volume fraction in Fig. [3 Note that this 
assumption was confirmed directly by a numerical exper- 
iment for a constant pressure system whose size is close 
to that of the system with the Lees-Edwards boundary 
conditions. Based on this assumption, we determine the 
volume fraction as a function of T at a constant pressure. 
In addition, using this, we can obtain the dependence 
of the shear stress on T at a constant pressure. Hence, 
we obtain (j/p as a function of I under constant n. For 
reference, we plot the lines I = const, in Fig [3] As an ex- 
ample in Fig. |3 let us consider the case that II = 1CP 3 
in which the volume fraction is larger than the critical 
one for a sufficiently small T. This case corresponds to 
the regime J < 10~ 3 and v > 0.650. Then, from Fig. 
121 we find that the shear stress remains almost constant. 
Next, in the interval 10~ 3 < J < 10 _1 , the states with 
II = 10 -3 are close to the critical line. Thus, in this 
regime, it is expected that a/p behaves as that in the 
critical line, and the scaling behavior given in Eq. is 
observed approximately. 

Generalizing the above discussion, we expect the typi- 
cal dependence of a/p on I as follows: 

a/p « const, for I <C J , (18) 

w I 1 / 5 for J < I <h, (19) 

when spatially homogeneous shear flow is realized. Note 
that Jo and I\ are dependent on the pressure. For ex- 
ample, for states with extreme pressures, such relations 
would not be observed. Since we wish to know the ex- 
tent to which this approximate power-law relation holds, 
in Fig. ^ we include the constant pressure data (J, a/p) 
obtained from Figs. [21 and |21 As expected from the 
above consideration, the power-law behavior is observed 
for the case in which II = 1.25 x 10 -3 . Furthermore, 
the system obeys the power-law regime even for the case 
IT = 5.4 x 10 -5 in which the line is located below the crit- 
ical states in Fig. |3 We do not understand the reason 
why the power law regime is so wide. 



V. CONCLUDING REMARK 

We have presented numerically and theoretically the 
scaling relations given in Eqs. J3J and Q for the system 



with the Lees-Edwards boundary conditions. From these 
new scaling relations, we also have an interpretation of 
the result observed in the system under constant pres- 
sure. The result is summarized in Eqs. I|18fl an d (|19fl . 

As far as we know, few experimental results exist in 
this regard. We expect that the power-law behaviors 
given in Eqs. J2Jl and Q are observed in systems with 
constant volume, e.g., by operating a rotating Coutte- 
flow system 0|. With regard to Eq. ((TBjl that is valid 
in the very low shear rate regime, we conjecture that a 
thermal activation process, which might lead to the log- 
arithmic dependence of a/p on J, occurs in this regime. 
Note that such behavior is ubiquitous in shear flow and 
sliding friction [l7|. Furthermore, the result reported in 
Ref. [HI might be related to Eq. JTHJ). We hope that 
more intensive experimental studies will be performed in 
this regard. 
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FIG. 5: a/p as a function of I for the model described in 
the caption of Fig. Qwith a Herzian type interaction fc'(l — 
Sr/(ri +r2)) 3//2 . The parameter values are as follows: L x /a — 
L y /a = 25, N = 10800, and r/^a/VkWn = 1.0. pa 2 /k' = 
3.75 x 10" 5 (data 1), 1.92 x 10" 3 (data 2). Furthermore, the 
square symbols represent the constant pressure data obtained 
from Figs. UjJandQ where the interaction force is k"((ri + 
r 2 ) - Sr/f /2 and p^a/k" = 6.625 x 10" 4 . 



Related to experimental studies, one may be interested 
in the dependence of our result on the choice of the model 
we investigate. For example, one may choose the Herzian 
type as an alternative for the interaction force. Indeed, 
such a model dependence has been discussed in the case 
of zero-temperature and zero applied stress 0. In our 
problem, it is highly expected that there is a critical state 
at which rheological properties exhibit power-laws. How- 
ever, it is not evident that the exponents remain the same 
values for the model with a Herzian type interaction. 

In order to consider the model dependence explicitly, 
we demonstrate the result of numerical experiments in 
Figs. H3 and [7| It is seen that the exponent of the 
system under constant pressure does not deviate so much 
from 1/5, while the exponents a and (3 seem to change 
slightly. We do not have a theoretical understanding for 
these values yet, because the time scale related to the 
particle collision is not directly determined from model 
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FIG. 6: Dimensionless shear stress E' as a function of 
the dimensionless shear rate F' for several values of v with 
the interaction force k'\(r\ + 7-2) — 5r) s ^ 2 . The dimension- 
less shear stress and the dimensionless shear rate are de- 
fined as £' = a^/a/k" and F' = 7^/ m/(fc"y / a), respec- 
tively. The parameter values are as follows: N = 1728, and 
t)/V^Wn = 1.0. 



parameters. It might be important to develop a theory 
by using a more physical time scale such as the collision 
interval. 

Finally, in our theoretical argument, Eq. Ijlll) plays 
an essential role. As in the case of dense colloidal sus- 
pensions |15|. one may investigate the pair-distribution 
function of granular systems in order to derive the order 
parameter equation. The establishment of a complete 



theory in which the scaling relations are derived from a 
microscopic model is an important topic for future stud- 
ies. 
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FIG. 7: Dimensionless shear stress II' as a function of the 
dimensionless shear rate V for several values of v with the 
interaction force k" ((n+r2) — Sr) s ^ 2 . The dimensionless shear 
stress and the dimen sionless she ar rate are defined as II' = 
p^/a/k" and T' = yy< m/(fc"y / a), respectively. The parameter 
values are as follows: N = 1728, and ri/^Jk" m,y/a = 1.0. 
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